close all;
clear all;

Cellsize = 5;       
fcount=10;
fig2=figure(2);

for i=1:fcount
    fname=fullfile('.',num2str(i),'/data_Cdc42T_time_course')
    clear bemt;
    bemt=importdata(fname,' ');
    if i==1
        Nf=size(bemt,1)/size(bemt,2);
        mcoord=zeros(fcount,Nf,2);
        N=size(bemt,2);
        dx  = Cellsize*sqrt(pi)/N; 
    end
    if mod(size(bemt,1),size(bemt,2))~=0
        display('Some frames were not fully written! The further progression of the script will be erroneous\n');
    else
        fig2=figure(2);
        for j=0:(Nf-1)
            bem=bemt(1+j*N:N+j*N,:);
            [xind,yind]=find(bem>mean(mean(bem)));
            mcoord(i,j+1,:)=[mean(xind)*dx, mean(yind)*dx];
        end
    end
end

figure(3);
hold on;
CM = jet(fcount);
for i=1:fcount
    scatter(mcoord(i,:,1),mcoord(i,:,2), 25,CM(i,:));
    plot(mcoord(i,:,1),mcoord(i,:,2),'color',CM(i,:));
    xlim([0 N*dx]);
    ylim([0 N*dx]);
    title('Movement of the geometric center of Cdc42T patch');
end
save mcoord_10runs_2.mat mcoord;
Nf=size(mcoord,2);
Nr=size(mcoord,1);
tf=30; %(min)
dt=tf/Nf;
msd=[];t=[];mstd=[];
for j=1:Nf
    sd_j=[];
    for k=1:Nr
        i=1;
        while (i+j)<=Nf
            dx=mcoord(k,i+j,1)-mcoord(k,i,1);
            dy=mcoord(k,i+j,2)-mcoord(k,i,2);
            sd_j=[sd_j (dx^2 + dy^2)];
            i=i+1;
        end
    end
    %hist(sd_j,0:0.01:max(sd_j));
    %size(sd_j)
    msd=[msd mean(sd_j)];
    mstd=[mstd std(sd_j)];
    t=[t dt*j/60];
end
figure(110);
plot(t,msd);
hold on;
errorbar(t,msd,mstd);
xlabel('Time, (hr)','fontsize',20);
ylabel('$$\langle \Delta X^2 \rangle$$, (micron)','interpreter','latex','fontsize',20);
title(['MSD for Bem geom center for ' num2str(tf) ' hr simulation'],'fontsize',26);
        